Table Of Contents

Previous topic

Programming examples I

Next topic

Testing Procedures

This Page

Programming examples II

250K Particles Fountain

Bellow I describe how to build the 250‘000 particles fountain demo included as example in the main application.

The demo can be executed with the command:

pyparticles_app -d fountain

First, you must import the libraries: numpy and particles set.

import numpy as np
import pyparticles.pset.particles_set as ps

For this demo I’ve chosen the midpoint method as ode numerical integration, but you can try with another one.

import pyparticles.ode.midpoint_solver as mds

The force is a combination between the constant gravity and the drag of the atmosphere, so I import the ‘const_force’ and ‘drag’; and the multiple force for building the full scenario

import pyparticles.forces.const_force as cf
import pyparticles.forces.drag as dr
import pyparticles.forces.multiple_force as mf

The ‘animation’ that control the simulation and the visualization is animation_ogl.

import pyparticles.animation.animated_ogl as aogl

And I use the default_boundary model to control the jet.

import pyparticles.pset.default_boundary as db

Let’s start

As first I’ve defined the function for modeling the jet. This function will be used by the boundary controller.

This is a time-dependent function. And for convenience time is defined as variable internal to the function. As rule, the time must be based on a reference to an object of type SimTime.

def default_pos( pset , indx ):

    t = default_pos.sim_time.time

    pset.X[indx,:] = 0.01 * np.random.rand( len(indx) , pset.dim )

    fs = 1.0 / ( 1.0 + np.exp( -( t - 2.0 ) ) )

    alpha = 2.0 * np.pi * np.random.rand( len(indx) )

    vel_x = 2.0 * fs * np.cos( alpha )
    vel_y = 2.0 * fs * np.sin( alpha )

    pset.V[indx,0] = vel_x
    pset.V[indx,1] = vel_y
    pset.V[indx,2] = 10.0 * fs + 1.0 * fs * ( np.random.rand( len(indx)) )

the main function

And finally you can start with the main function:
As first I’ve defined the number of steps the dt (step time), the size of the particle set and the particles set pset
def fountain():

    steps = 10000000
    dt = 0.01

    pcnt = 250000

    pset = ps.ParticlesSet( pcnt )

Initial position

As a second point define the initial positions of the particles, in order to generate a continuous stream.

pset.M[:] = 0.1
pset.X[:,2] = 0.7 * np.random.rand( pset.size )

Forces model

The force model is based on a combination between a constant force, for simulating the gravity, and the drag for simulating the friction of the air.

The two forces must be stored in an object of type MultipleForce.

grav = cf.ConstForce( pset.size , dim=pset.dim , u_force=( 0.0 , 0.0 , -10.0 ) )
drag = dr.Drag( pset.size , dim=pset.dim , Consts=0.01 )

multi = mf.MultipleForce( pset.size , dim=pset.dim )

multi.append_force( grav )
multi.append_force( drag )

multi.set_masses( pset.M )

ODE integration

As a ODE numerical integration method I’ve used the mid point algorithm.

solver = mds.MidpointSolver( multi , pset , dt )
solver.update_force()

Simulation time

Set up the simulation time in the default_pos function, used for modeling the jet.

default_pos.sim_time = solver.get_sim_time()

Build the boundary model

The tuples bd represents the size of the box closed domain: ( min_x , max_x , min_y , max_y , min_z , max_z )

Note

The class DefaultBoundary positions the particles exited from the limits of the domain according to the function defualt_pos

bd = ( -100.0 , 100.0 , -100.0 , 100.0 , 0.0 , 100.0 )
bound = db.DefaultBoundary( bd , dim=3 , defualt_pos=default_pos )

pset.set_boundary( bound )

Build the ‘animation’ class and start

Note

a.init_rotation( -80 , [ 0.7 , 0.05 , 0 ] )
setup an initial rotation where parameters are ( rot angle , axis of rotation )
a.draw_particles.set_draw_model( 1 )
enable the vectorized rendering, that is fundamental for drawing 250‘000 particles
a = aogl.AnimatedGl()

a.ode_solver = solver
a.pset = pset
a.steps = steps

a.draw_particles.set_draw_model( 1 )

a.init_rotation( -80 , [ 0.7 , 0.05 , 0 ]  )

a.build_animation()
a.start()

return